In [56]:
import numpy as np
import pylab as plt
import isqg #the isQG package
#define longitude and latitude for the regional domain
lon = np.linspace(0,10,50)
lat = np.linspace(30,40,50)
#transform the coordinates from degrees to meters
x,y = isqg.lonlat2xy(lon,lat)
#define a gaussian function
gauss = lambda x,y,r: np.exp(- (
( x-x.mean() )**2 + ( y-y.mean() )**2
) /2./r**2 )
#set the surface density perturbation
radi = 1e5 # the radius of the gaussian eddy
ssd = 0.1*gauss(x,y,radi) #surface density anomaly of the gaussian eddy
#vertical levels and stratification profile N^2
z = np.linspace(0,-2000,30)
n2 = np.ones_like(z)* 1e-5
# technically you can set ssh here, but for a demonstration purpurse
# we will set SSH based on $\psi^s$ which will be calculated later
ssdm = np.ones_like(ssd)*ssd.mean()
ssda = ssd - ssd.mean()
# reconstruct mean density based on mean surface density and N^2
rhom = isqg.N2rho(ssdm,isqg.twopave(n2), np.diff(z))
# initialize isqg data
d = isqg.isqg_data(n2=n2,z=z,rhosorg=ssda,lon=lon,lat=lat,rho0=1035.)
# set the bottom boundary for SQG solution, default is 'rho=0'.
# 'psi=0' is no bottom velocity condition
d.bottomboundary='psi=0'
d.solve_sqg() # SQG solution is calcuated after this call
In [57]:
plt.figure(figsize=(10,10))
plt.subplots_adjust(wspace=0.2,hspace=0.2)
plt.subplot(221)
plt.contourf(ssda,levels=np.arange(-0.06,0.11,0.02))
plt.colorbar()
plt.title('$\\rho^o$',fontsize=30)
plt.subplot(222)
plt.contourf(d.rhos[0,:,:],levels=np.arange(-0.06,0.11,0.02))
plt.colorbar()
plt.title('$\\rho^s$',fontsize=30)
plt.subplot(223)
plt.contourf((ssda-d.rhos[0,:,:])/ssda.max()*100)
plt.colorbar()
plt.title('$\\rho^o-\\rho^s$',fontsize=30)
plt.subplot(224)
plt.contourf(lon,z,d.rhos[:,:,25])
plt.colorbar()
plt.title('a vertical section')
show()
In [58]:
# add surface height field
# the ssh is similar to the density induced ssh but has 1 degree eastward shift
d.ssh = np.roll(1.2*d.psis[0,:,:]*d.f0/9.81,5,axis=1)
#this call solves interior solution
d.solve_psii()
In [59]:
levs = np.arange(-8000,26000, 2500)
plt.figure(figsize=(10,10))
plt.subplots_adjust(wspace=0.2,hspace=0.2)
plt.subplot(221)
plt.contourf(d.psit[:,:], levels=levs)
plt.colorbar()
plt.title('$\psi^t$',fontsize=30)
plt.contour(d.psit, levels=levs, colors='w')
plt.subplot(222)
plt.contourf(d.psis[0,:,:], levels=levs)
plt.colorbar()
plt.title('contour:$\psi^t$;color:$\psi^s$',fontsize=20)
plt.contour(d.psit, levels=levs, colors='w')
plt.subplot(223)
plt.contourf(d.psii[0,:,:], levels=levs)
plt.colorbar()
plt.title('contour:$\psi^t$;color:$\psi^i$',fontsize=20)
plt.contour(d.psit, levels=levs, colors='w')
plt.subplot(224)
plt.contourf(d.psii[0,:,:]+d.psis[0,:,:], levels=levs)
plt.colorbar()
plt.title('contour:$\psi^t$;color:$\psi^i +\psi^s$',fontsize=20)
plt.contour(d.psit, levels=levs, colors='w')
plt.show()
In [60]:
#use dir(d) to check other available variables
dir(d)
Out[60]: